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Semiconductor  Measurement  Technology: 
Version  2.0  of  the  TXYZ  Thermal  Analysis  Program:  TXYZ20 

John  Albers 
Semiconductor  Electronics  Division 
National  Institute  of  Standards  and  Technology 
Gaithersburg,  Maryland  20899 

Abstract 

The  TXYZ  computer  program  has  been  used  for  a  number  of  years  for  the  thermal  analysis 
of  semiconductor  devices  and  packages.  This  program  makes  use  of  the  closed  form,  Fourier 
series  solution  of  the  steady-state  heat  flow  equation  for  the  general  case  of  a  rectangular 
three-layer  structure  with  multiple  heat  sources  on  the  top  surface.  TXYZ  provides  for  the 
calculation  of  the  temperature  at  any  set  of  points  in  this  structure  and  has  proven  useful 
for  the  determination  of  the  steady-state  temperature  distribution  of  semiconductor  chips 
and  packages.  This  report  presents  TXYZ  Version  2.0  (TXYZ20)  which  is  a  revised  and 
updated  version  of  the  original  TXYZ  program.  The  TXYZ20  program  incorporates  more 
flexible  handling  of  input  data,  assignment  of  positive  or  negative  noninteger  weights  to 
the  various  heat  sources  or  heat  sinks,  and  improved  evaluation  of  limiting  forms  in  the 
code. 

The  first  part  of  this  report  consists  of  a  discussion  of  the  general  elements  in  the  TXYZ 
code  and  the  particular  changes  which  have  been  made  to  it  to  obtain  TXYZ20.  The 
second  part  of  the  report  contains  a  discussion  of  several  examples  of  the  running  of  the 
code.  Several  annotated  input  data  files  are  presented  and  discussed  to  show  both  the 
increased  flexibility  of  the  input  data  and  the  actual  use  of  the  updated  code.  Running  the 
TXYZ20  code  for  one  of  the  input  files  provides  a  benchmark  for  several  machines.  The 
user  may  wish  to  run  this  example  for  the  purpose  of  comparing  the  CPU  times  involved. 
The  appendix  contains  an  annotated,  internally  documented  listing  of  the  FORTRAN 
source  code  for  TXYZ20. 

The  FORTRAN  source  code  (total  of  about  21  kbytes)  and  sample  input  and  output  data 
files  are  available  in  ASCII  format  using  a  number  of  transfer  vehicles.  These  include: 
standard  8-track  magnetic  tape  (ASCII,  density  =  1600,  record  =  80,  block  =  1600),  5.25- 
in.  (360-kbyte  and  1.2-Mbyte)  DOS-formatted  floppy  disks,  and  electronic  mail  over  the 
Internet.  The  sample  input  and  output  data  files  are  included  so  that  the  user  can  check 
the  program  for  proper  operation  as  well  as  to  become  acquainted  with  the  setup  and  use 
of  the  code.  Users  of  the  TXYZ  code  will  find  the  updated  TXYZ20  code  easy  to  use  and 
should  benefit  from  the  more  flexible  input  and  the  more  general  treatment  of  heat  sources 
and  heat  sinks. 

Key  words:  FORTRAN;  Fourier  analysis;  integrated  circuit;  semiconductor  devices;  semi- 
conductor materials;  steady-state  heat  flow;  thermal  analysis;  thermal  conductivity. 
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Introduction 


The  operation  of  a  semiconductor  device  relies  on  the  passage  of  current  through  portions 
of  the  structure.  This  is  accompanied  by  power  dissipation  and  heating  which  gives  rise 
to  a  temperature  distribution.  In  time,  thermal  stresses  may  arise  and  lead  to  device 
degradation  and  possible  failure.  In  addition  to  the  effects  of  thermal  stresses  on  device 
reliability,  the  modeling  of  the  steady-state  thermal  response  of  a  system  is  also  useful 
since  the  temperature  distribution  may  have  an  effect  upon  the  local  mechanical,  thermal, 
electrical,  or  optical  properties  of  the  system.  Consequently,  an  accurate  physical  model 
of  the  temperature  distribution  under  the  power  condition  of  actual  operation  is  of  great 
importance.  It  is  also  possible  to  ascertain  the  relative  effects  of  composition  (thermal 
conductivity)  and  geometry  (layer  thickness).  By  carefully  optimizing  composition  and 
geometry,  it  may  be  possible  to  minimize  the  thermal  stress  and  hence  ensure  optimal 
device  lifetime. 

The  physical  and  mathematical  model  used  here  is  taken  from  the  work  of  Kokkas  [1].  This 
has  been  used  previously  to  construct  the  TXYZ  program.  The  Kokkas  model  and  the 
TXYZ  code  have  been  described  in  considerable  detail  in  NBS  Special  Publication  400-76 
[2].  This  report  is  available  and  should  be  consulted  for  the  details  of  the  analysis  which 
are  still  current  and  valid. 

Review  of  the  Model,  the  TXYZ  Code,  and  the  TXYZ20  Code 

For  the  purposes  of  providing  a  background  for  the  TXYZ20  code,  a  brief  outline  of  the 
salient  features  of  the  earlier  report  on  TXYZ  is  presented  here.  It  is  assumed  that  the 
reader  has  a  copy  of  reference  [2].  The  symbols  used  here  are  denned  in  this  reference. 
The  outline  presented  here  should  provide  the  reader  with  an  overview  of  the  analysis  and 
a  quick  guide  as  to  where  to  find  the  various  portions  of  the  calculation. 

The  starting  point  of  the  calculation  is  the  steady-state  heat  flow  equation  for  a  single 
rectangular  layer.  In  Cartesian  coordinates,  this  takes  on  the  form 

r    Q2  Q2  Q2  ^| 

This  equation  is  solved  in  this  coordinate  system  using  Fourier  analysis.  The  requirement 
that  there  be  no  heat  transport  out  of  the  lateral  boundaries  of  the  structure  leaves  only 
"quantized"  cosine  functions  in  the  Fourier  representation  of  the  temperature, 

V^1        4r(n,m,  z)  cos(mvx/ Lx)  cos(rmry / Ly) 
T{x,V,z)  =  ^  2-,  (*n0  +  i)(*m0  +  1)LxL  ' 

n=0  m=0  v  ' x  '  9 

where  8nn>  is  the  Kronecker  delta  and  is  equal  to  unity  if  n  =  n'  and  zero  otherwise. 
The  Fourier  coefficients,  r(n,m,z),  are  given  by 

r(n,  m,  z)  =  a  cosh(72)  +  j3  sinh(7-z), 
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where  the  coefficients  a  and  /?,  which  may  be  functions  of  7,  are  determined  from  the  two 
z-dependent  boundary  conditions. 


The  above  analysis  provides  the  basis  for  the  general  multiple-layer  problem  in  terms  of 
the  Fourier  coefficients.  For  this  case,  it  is  assumed  that  the  temperature  satisfies  the 
steady-state  equation  in  each  of  the  layers.  Consequently,  in  each  of  the  layers,  the  Fourier 
expansion  of  the  temperature  takes  the  form 

T{(n,  m,  z)  =  oti  cosh(7z)  +     sinh(7z),  (i  =  1,  ...,n). 

For  each  layer,  there  are  two  coefficients  which  arise  from  the  fact  that  the  basic  tem- 
perature equation  is  second  order.  These  are  determined  from  the  z-dependent  boundary 
conditions.  The  general  set  of  z-dependent  boundary  conditions  for  the  multiple  layer 
problem  is  provided  by  the  conditions  that:  heat  enters  (or  leaves)  the  structure  through 
the  portions  of  the  top  layer  where  power  is  applied,  the  temperature  and  the  heat  flow 
are  continuous  across  the  internal  interfaces,  and  the  temperature  is  continuous  across  the 
interface  between  the  bottom  layer  and  the  heat  sink.  The  above  equations  and  boundary 
conditions  provide  the  mathematical  framework  in  which  the  multiple  layer  problem  may 
be  solved.  For  the  purposes  of  numerical  modeling,  the  multiple  layer  problem  is  special- 
ized to  one  with  three  layers.  In  this  case,  the  Kokkas  [1]  analysis  leads  to  the  solution  in 
the  form 

,         \      p  4U(n,Tn)Tj(n,m,  z)  cos(nirx/Lx)  cos(mny/Ly) 

lt{x,y,z)  _  f0  2.  (*„o  +  l)(*mo ' 


n=0  m=0 


The  symbols,  analysis  and  particular  forms  of  Ti(n,m,  z)  are  discussed  in  some  detail 
in  reference  [2].  There  are  several  parts  of  the  above  equation  which  can  be  identified. 
First,  there  is  the  z-dependent  part  of  the  solution  in  each  layer,  Ti(n,m,  z).  The  two 
cosine  functions,  cos(mvx  /  Lx)  cos(rmvy / Ly),  arise  from  the  x-  and  y-dependent  parts  of 
the  steady-state  equation.  Their  specific  form  arises  from  the  requirement  of  the  lateral 
boundary  conditions.  Finally,  the  requirement  that  heat  enters  or  leaves  the  top  surface 
only  through  the  heating  elements  introduces  the  function  PoU(n,m).  It  should  be  noted 
that  Pq  appears  only  in  front  of  the  sum.  It  is  a  scale  or  multiplicative  factor  which  for 
convenience  may  be  set  equal  to  unity.  Its  effect  is  to  scale  the  temperature  in  a  direct 
linear  fashion.  Performing  the  calculation  for  PQ  =  1  for  a  given  structure  then  provides 
the  temperature  distribution  for  arbitrary  P0  which  can  then  be  obtained  by  multiplying 
by  the  Po  which  is  used. 

In  the  above,  the  function,  ?7(n,m),  is  related  to  the  heat  sources  (sinks)  on  the  top  surface 
of  the  structure.  This  function  is  defined  as  the  double  cosine  transform  of  the  geometric 
part  of  the  surface  power  function,  U(x,y), 

U(n,m)  =  /       /  U(x,y)cos(mvx/Lx)cos(mny/Ly)dxdy. 
Jo  Jo 

The  detailed  form  of  this  function  is  considered  in  reference  [2].  There,  the  relation  of 
this  function  to  the  location  and  size  of  the  heat  sources  was  presented.  In  addition,  the 


limiting  form  of  the  function  as  either  n  or  m  approach  zero  was  considered,  and  the  correct 
limit  was  established  for  the  purposes  of  numerical  analysis.  The  reader  is  referred  to  that 
report  for  details. 

In  the  original  version  of  TXYZ,  all  heat  sources  had  equal  unit  weights.  If  the  user  wanted 
to  construct  a  more  general  heat  source,  this  could  only  be  done  in  unit  increments  by 
superimposing  sources.  Hence,  to  perform  the  calculation  for  a  heating  element  which  had 
five  times  the  effect  as  another,  the  element  would  have  to  constructed  five  times.  Only 
integer  increments  were  possible.  The  first  change  made  in  TXYZ  is  the  generalization 
of  the  heating  elements.  It  should  be  noted  that  this  does  change  the  basic  model.  The 
updated  TXYZ20  program  allows  for  noninteger  weights  to  be  assigned  to  each  heating 
element.  This  weight  may  be  positive  or  negative  and  hence  represent  the  effects  of  either 
a  heat  source  or  a  heat  sink.  The  weight  is  entered  through  the  input  data  file  and  appears 
with  the  location  and  size  of  the  element.  The  reader  should  refer  to  the  last  lines  in 
the  input  data  files  presented  in  tables  1,  4,  5,  and  7  for  the  use  of  weighting  of  the  heat 
elements.  The  more  general  information  is  then  passed  from  the  main  program  to  the 
UZERO  function  through  a  COMMON  statement.  The  weighting  of  the  heating  element 
was  used  originally  in  a  calculation  of  the  thermal  interactions  between  electromigration 
test  structures  [3].  The  careful  and  consistent  evaluation  of  the  temperature  is  important  in 
this  case  as  the  test  structures  are  used  in  accelerated  testing  and  the  temperature  profile 
generated  by  one  of  the  elements  may  increase  the  temperature  at  adjoining  elements 
which  may  accelerate  the  failure  rates.  As  electromigration  is  temperature  dependent,  it 
is  of  central  importance  to  evaluate  the  temperature.  This  is  especially  the  case  where  the 
thermal  interaction  of  neighboring  heating  elements  comes  into  play. 

The  weight  is  entered  on  the  same  line  of  the  input  data  file  where  the  position  and  lengths 
along  x  and  y  are  entered.  Instead  of  having  to  enter  the  same  element  five  times,  a  weight 
of  5  may  simply  be  used.  In  addition,  noninteger  weights  which  are  positive  or  negative 
may  be  entered  to  represent  sources  or  sinks.  This  feature  should  increase  the  versatility 
of  the  code  to  a  more  general  class  of  problems. 

Returning  to  the  overview  of  the  original  TXYZ  code  and  the  basis  for  the  evolution  of 
TXYZ20,  the  behavior  of  the  Fourier  coefficients  for  both  small  and  large  values  of  the 
arguments  is  next  on  the  agenda.  The  importance  of  this  aspect  of  the  analysis  and 
the  resulting  numerical  code  is  found  in  the  fact  that  the  equations  involved  contain  the 
possibility  of  numerical  infinities  or  overflow  situations.  These  occur  for  n  =  0  and  m  =  0 
as  well  as  for  large  values  of  these  summation  variables.  The  reason  for  the  relatively 
compact  nature  of  the  original  TXYZ  code  and  its  numerical  efficiency  is  that  great  care 
was  exercised  in  the  investigation  of  the  functions  for  the  small  and  large  argument  regime. 
For  small  values  of  the  argument,  the  evaluation  of  limiting  forms  led  to  expressions  which 
were  correct  and  free  of  artificial  numerical  singularities.  For  large  values  of  the  argument, 
special  care  was  required.  Many  of  the  functions  involved  contained  the  hyperbolic  sinh  and 
cosh  functions.  These  approach  the  exponential  function  for  large  values  of  the  argument. 
These  provide  the  possibility  for  numerical  infinities  or  overflows.  Careful  investigation 
of  the  forms  of  the  functions  showed  that  these  numerical  infinities  were  removable  by 
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other  numerical  infinities.  For  example,  an  exponential  in  the  numerator  and  denominator 
of  a  function  would  give  rise  to  two  overflows  when  evaluated  by  a  computer  but  would 
cancel  each  other  and  give  rise  to  a  finite  result  when  analytic  evaluation  was  properly 
applied.  The  construction  of  the  function,  FUNZ,  in  the  TXYZ  code  was  based  upon  such 
evaluations  and  their  numerical  implementations.  With  the  passage  of  time  and  the  use 
of  the  TXYZ  code,  one  place  where  such  cancellation  was  found  to  be  incomplete  was  in 
the  case  of  a  fairly  thick  middle  layer.  For  most  cases  of  interest  to  the  modeling  of  the 
steady-state  thermal  response  of  semiconductor  device  structures,  these  layers  are  usually 
thin.  However,  this  situation  has  now  been  carefully  investigated  and  the  numerically 
stable  limiting  form  has  been  determined  and  used  in  the  code  in  the  middle  layer  portion 
of  the  FUNZ  function.  The  introduction  of  the  nonsingular  limiting  form  is  contained  in 
the  TXYZ20  code. 

The  original  TXYZ  report  contains  a  section  showing  how  the  one- dimensional  equations 
are  obtained  in  the  special  case  of  a  power  source  which  completely  covers  the  top  surface. 
This  analysis  requires  the  careful  evaluation  of  limits.  When  properly  applied,  the  one- 
dimensional  expressions  for  the  thermal  resistance  are  obtained.  The  reader  is  referred  to 
reference  [2]  for  the  details. 

The  last  item  to  discuss  before  addressing  the  structure  of  the  input  data  file  for  TXYZ20 
is  concerned  with  the  effect  of  the  upper  summation  limits  on  the  calculated  temperature. 
While  this  is  discussed  on  pp.  32-33  and  in  figure  7  of  reference  [2],  it  is  an  important 
concept  to  understand  and  utilize  correctly  when  running  the  code.  Consider  the  case 
of  a  structure  with  lateral  dimension  of  L  and  a  heat  element  of  lateral  dimension  of  A. 
For  the  Fourier  series  to  begin  "seeing"  this  element,  the  cosine  function  in  the  basis  set 
must  have  at  least  one  complete  cycle  in  A.  As  the  element  is  rectangular,  one  cycle  is 
certainly  not  enough,  as  the  cosine  is  a  poor  representation  of  the  rectangle.  Consequently, 
higher  "frequencies"  are  in  order  to  assure  the  adequacy  of  the  representation.  The  rule 
of  thumb  is  that  the  number  of  terms  should  be  at  least  on  the  order  of  L/ A.  A  stronger 
rule  of  thumb  would  require  several  times  the  above  rule  of  thumb.  Certainly,  the  example 
of  reference  [2]  for  L/A  =  200  shows  that  the  asymptotic  value  of  the  temperature  is 
attained  for  400  to  500  terms.  In  general,  the  stronger  rule  of  thumb  should  be  applied  to 
the  smallest  heat  element  in  order  to  achieve  better  accuracy. 

In  its  present  form,  TXYZ20  allows  an  upper  limit  of  500  terms.  With  smaller  feature 
sizes,  this  upper  limit  may  be  inadequate.  If  this  is  the  case,  it  is  important  to  edit  the 
FORTRAN  source  code  for  TXYZ20.  This  requires  the  following  changes.  First,  the 
DIMENSION  statements  in  the  main  program  for  COSTY(500),  ARUZER(500,500)  and 
ARFUNZ(500,500)  must  be  increased.  Second,  the  line  which  tests  the  input  values  of  the 
two  upper  limits  and  stops  execution  of  the  program  if  they  are  too  big  must  be  turned 
off.  The  evolution  of  computers  over  the  past  decade  has  put  more  memory  in  the  hands 
of  the  user.  Increasing  the  size  of  the  DIMENSION  statements  may  just  test  the  extent 
of  this  memory.  Failure  of  the  modified  code  to  compile  due  to  memory  restrictions  is  the 
test. 

A  note  of  caution  is  in  order  here.  If  the  test  line  is  turned  off  and  the  upper  limit  is  larger 
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than  the  size  of  the  DIMENSION  statement,  the  calculation  will  proceed  with  values  being 
stored  in  incorrect  registers.  No  meaning  should  be  ascribed  to  the  results.  Of  course,  the 
author  is  available  for  discussion  and  advice. 

Examples  of  Calculations  with  TXYZ20 

This  portion  of  the  report  is  concerned  with  examples  of  the  input  data  file  for  the  TXYZ20 
code.  The  reader  is  encouraged  to  read  through  the  listing  of  the  FORTRAN  source  code 
for  TXYZ20.  In  regard  to  the  input  data,  special  attention  should  be  given  to  the  internal 
documentation  which  is  provided  in  the  portion  which  describes  the  input  section  (see  pp. 
21-22).  The  variables  and  their  use  as  well  as  the  form  of  their  input  are  all  described  in 
considerable  detail  there.  Before  describing  the  changes  in  the  input  data  files,  a  special 
note  is  in  order  concerning  the  units  used  in  the  calculations.  It  is  of  utmost  importance 
that  the  various  input  parameters  be  in  a  consistent  set  of  units.  Whatever  length  unit  is 
used  for  the  dimensions  of  the  structure,  the  dimensions  of  the  heat  elements  and  the  layer 
thicknesses  must  also  be  used  in  the  thermal  conductivities  and  the  power  density.  Hence, 
if  length  is  used  for  the  former,  then  the  thermal  conductivities  must  be  in  the  units  of 
watts/  °C  length  and  the  power  density  must  be  in  the  units  of  watts/(length)2. 

The  major  changes  in  the  input  have  to  do  with  (1)  the  input  or  generation  of  the  x,  y, 
and  z  values  to  be  used  and  (2)  the  input  of  the  weighting  factor  for  the  heating  elements. 
These  are  discussed  in  some  detail  here. 

The  input  for  the  x,  y,  and  z  values  has  been  made  more  flexible  and  consistent.  As  the 
treatment  of  each  of  these  is  the  same,  discussion  will  be  focused  on  only  the  first. 

The  first  five  lines  of  the  input  data  file  remains  the  same  as  that  used  in  the  original 
TXYZ  code.  These  contain  the  dimensions  of  the  rectangular  structure,  the  thickness,  and 
thermal  conductivity  of  the  three  layers  and  the  summation  limits  for  the  double  sum.  The 
sixth  line  contains  the  integer  index,  IEDGEX.  This  is  used  for  the  generation  or  input  of 
the  values  of  x  to  be  used  in  the  calculation.  If  IEDGEX  =  1,  then  a  set  of  equispaced 
values  of  x  will  be  generated.  For  this  case,  the  next  line  of  the  input  data  file  is  to  contain 
three  numbers.  These  are:  ILX  (the  number  of  x  values  to  be  generated),  XI  (the  value 
of  the  first  x  point),  and  STEPX  (the  value  of  the  increment  in  x).  Once  these  numbers 
are  read  in,  the  program  constructs  the  set  of  points  according  to 

X(I)  =  XI  +  (I-1)*STEPX  (1=1,..., ILX). 

If  IEDGEX  =  2,  then  a  set  of  values  of  x  will  be  read  from  the  data  file.  For  this  case,  the 
the  next  line  of  the  input  data  file  will  contain  the  number  of  x  values  to  be  read,  ILX. 
The  program  will  then  read  in  ILX  lines  of  input  data  file  containing  one  value  of  x  for 
each  line.  The  IEDGEX  =  2  option  allows  for  the  calculation  to  be  performed  at  a  set  of 
x  values  which  are  not  necessarily  equispaced.  This  is  particularly  useful  if  the  user  wants 
to  have  some  of  the  x  values  closely  spaced  on  or  near  a  heat  element  and  other  values 
spread  out  away  from  the  heat  element.  This  option  is  motivated  by  the  fact  that  the 
temperature  is  usually  rapidly  varying  near  the  element  and  slowly  varying  away  from  it. 
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Once  the  values  of  x  have  been  generated  or  read,  the  code  goes  on  to  generate  or  read 
the  values  of  y  and  then  to  generate  or  read  the  values  of  z  by  the  same  construction.  As 
noted  in  the  source  code  internal  documentation  (see  pp.  21-22),  the  information  on  the 
generation  or  reading  of  the  z  values  is  always  entered  as  positive  quantities.  The  code 
takes  care  of  making  the  appropriate  sign  changes. 

The  other  major  change  in  the  input  data  file  is  the  introduction  of  the  weighting  factor 
for  each  heat  element.  Once  the  number  of  heat  elements  has  been  read  in  along  with 
the  power  density,  the  information  for  each  of  the  heating  elements  is  then  read.  For  each 
element,  a  single  line  of  five  numbers  contains  all  of  the  information.  The  weighting  factor 
is  the  first  entry  on  the  line.  Hence,  the  line  contains  WTSOUR(I),  XSOUR(I),  YSOUR(I), 
LXSOUR(I),  AND  LYSOUR(I).  WTSOUR(I)  is  the  weighting  factor  (positive  for  source, 
negative  for  sink)  for  the  i-th  element.  This  is  a  REAL  number  and  does  not  have  to  be 
INTEGER. 

The  above  features  make  the  structure  of  the  input  data  file  more  flexible  and  consistent. 
It  is  important  to  note,  as  has  been  discussed  before,  that  the  power  density  variable, 
Po,  be  set  equal  to  unity.  In  effect,  the  calculation  then  generates  a  scaled  temperature, 
T(x,y,x)/Po,  which  may  always  be  rescaled  for  a  particular  value  of  Po. 

Before  turning  to  the  examples,  it  is  also  important  to  note  the  changes  made  in  the 
output  data  file.  The  format  used  for  the  construction  of  the  output  data  is  of  a  general 
nature.  It  has  been  written  to  accommodate  a  wide  range  of  numerical  values  without  loss 
of  generality  or  data.  This  was  motivated  by  the  fact  that  datum  may  be  lost  on  output 
if  it  exceeds  the  format  of  the  write  statement.  In  order  to  circumvent  this,  the  data  may 
appear  to  have  more  significant  figures  than  might  be  reasonable  or  needed.  The  user 
should  keep  this  in  mind  when  looking  at  the  output.  The  user  should  round  the  output 
following  accepted  rules  for  this  process.  All  the  WRITE  statements  in  the  FORTRAN 
source  code  direct  output  to  logical  device  for012.  The  one  exception  is  the  line  which  is 
written  if  either  one  or  both  of  the  upper  summation  limits  exceeds  the  dimensionality  of 
the  arrays  used.  This  terminates  program  execution.  Assuming  that  this  is  not  the  case, 
the  data  which  are  written  to  for012  are  x,  y,  z,  T(x,y,z).  The  major  loop  in  the  code 
(3000)  is  involved  with  the  various  z  values  used.  Once  the  list  of  x,  y,  and  z  values  for  the 
calculations  of  T(x,y,z)  has  been  completed,  the  remaining  section  of  the  output  data  file 
contains  the  listing  of  the  input  variables  used  in  the  calculations.  Saving  the  output  data 
without  some  identifier  of  the  input  data  can  present  problems  as  to  where  the  output 
came  from.  This  is  taken  care  of  with  the  present  listing  to  for012. 

Benchmark  Example 

The  first  example  is  contained  in  the  two  files,  inputl.dat  and  outputl.dat.,  This  is  for  the 
case  of  a  single  heat  source  of  unit  weight.  The  temperature  is  calculated  for  11  values  of 
x.  The  example  is  illustrative  of  the  use  of  the  TXYZ20  code.  The  reader  is  referred  to 
the  annotated  input  and  output  data  files  (see  tables  1  and  2)  where  each  line  of  the  file 
is  described.  This  particular  problem  is  of  interest  as  it  requires  a  small  amount  of  time 
to  carry  out  but  does  provide  for  some  indication  of  the  relative  speeds  which  users  can 
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Table  1.    Annotated  input  data  file  (inputl.dat)  for  TXYZ20  code. 

The  data  correspond  to  those  contained  in  the  file  inputl.dat  which  can  be 
used  to  run  the  case  of  a  single  heat  source  of  unit  weight.    This  illustrates 
the  use  of  the  IEDGEX  =  1  generation  of  a  set  of  values  of  x  to  be  used  in  the 
calculation.     In  this  and  subsequent  annotated  input  data  file  listings,  it 
is  important  for  the  user  to  keep  in  mind  the  comments  of  p.  6  regarding  the 
consistency  of  the  units  to  be  used  in  the  calculation. 


1    40.0    20.0    40.0    20.0     WEIGHT,  X,  LENGTH  ALONG  X,  Y,  LENGTH  ALONG  Y 


INPUT  DATA 


DESCRIPTION  OF  INPUT  DATA 


100  100 
1  .001 
5  .0001 
4  .005 
500  500 
1 


11  0.0  10. 

1 

1  50.  0.0 
1 

1  0.0  0.0 

1  1.0 


X  AND  Y  DIMENSIONS  OF  RECTANGULAR  STRUCTURE 
THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  TOP  LAYER 
THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  MIDDLE  LAYER 
THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  BOTTOM  LAYER 
UPPER  SUMMATION  LIMITS  FOR  N  AND  M  SUMMATIONS 
IEDGEX  (=1,  THEN  READ  IN  THREE  VALUES  ON  NEXT  LINE) 
NUMBER  OF  VALUES,  FIRST  POINT,  INCREMENT 
IEDGEY  (=1,  THEN  READ  IN  THREE  VALUES  ON  NEXT  LINE) 
NUMBER  OF  VALUES,  FIRST  POINT,  INCREMENT 
IEDGEZ  (=1,  THEN  READ  IN  THREE  VALUES  ON  NEXT  LINE) 
NUMBER  OF  VALUES,  FIRST  POINT,  INCREMENT 
NUMBER  OF  SOURCES  AND  POWER  DENSITY 
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Table  2.    Annotated  output  data  file  (outputl.dat)  for  TXYZ20  program. 

The  data  correspond  to  those  contained  in  the  file  outputl.dat  which  is 
generated  for  the  case  of  a  single  heat  source  of  unit  weight. 
The  first  block  of  lines  are  the  x,y,z,  and  calculated  t(x,y,z)  values. 
Beyond  these  are  the  listing  of  the  input  variables.    In  this  and  subsequent 
annotated  output  data  file  listings,   it  is  important  for  the  user  to  keep  in 
mind  the  comments  of  p.  7  regarding  the  format  used  for  the  output  data  and 
the  possible  need  for  rounding. 


X 

y 

z 

t(x,y,z) 

.0000000E+00 

50 

00000 

0 

0000000E+00 

107.2166 

10 

00000 

50 

00000 

0 

0000000E+00 

229.7415 

20 

00000 

50 

00000 

0 

0000000E+00 

899.0788 

30 

00000 

50 

00000 

0 

0000000E+00 

3868.299 

40 

00000 

50 

00000 

0 

0000000E+00 

17908.64 

50 

00000 

50 

00000 

0 

0000000E+00 

29526.74 

60 

00000 

50 

00000 

0 

0000000E+00 

17908.62 

70 

00000 

50 

00000 

0 

0000000E+00 

3868.297 

80 

00000 

50 

00000 

0 

0000000E+00 

899.0831 

90 

00000 

50 

00000 

0 

0000000E+00 

229.7416 

100.0000 

50 

00000 

0 

0000000E+00 

107.2136 

STEADY-STATE  THERMAL  ANALYSIS  CALCULATION  USING  EQS.  (13) -(23)  OF  KOKKAS 

THERMAL  CONDUCTIVITIES  AND  LAYER  THICKNESSES 
LX=    100.00       LY=  100.00 

Ll=  1.00000  L2=  5.00000  L3=  4.00000 
Kl=  0.00100000  K2=  0.00010000  K3=  0.00500000 

UPPER  SUMMATION  LIMITS       NUP=    500    MUP=  500 


NUMBER  OF  HEAT  S0URCES=  1 
POWER  DENSITY=  1.000000 

COORDINATES,  LENGTHS,  AND  WIDTHS  OF  HEAT  SOURCES 

HEAT  SOURCE       WTSOUR        XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000         40.00000         40.00000         20.00000  20.00000 
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expect  to  see  in  the  execution  of  the  code.  This  particular  input  data  file  was  used  for  the 
TXYZ20  code  on  several  machines  [4].  These  include:  a  VAX  11/785,  a  SUN,  a  SPARC  2 
and  a  CRAY.  All  of  these  systems  had  sufficient  memory  to  compile  the  code  and  run  it 
for  the  upper  limits  of  500.  Users  of  the  code  on  PCs  may  have  to  experiment  to  see  if  the 
dimension  statements  and  the  upper  limits  are  too  large.  This  will  be  influenced  by  the 
amount  of  RAM  as  well  as  the  presence  of  a  coprocessor.  The  execution  of  the  code  for 
this  problem  shows  a  range  of  relative  speeds  from  1  to  55.  The  results  of  this  benchmark 
example  are  presented  in  table  3. 

The  above  example  makes  use  of  the  IEDGE=1  option  for  data  entry.  This  can  also 
be  accomplished  by  using  the  IEDGE=2  option.  This  alternative  is  contained  in  the  file 
inputla.dat  which  is  annotated  in  table  4.  The  output  from  running  these  data  is  the  same 
as  contained  in  outputl.dat  (illustrated  in  table  2). 

Other  Examples 

Several  other  examples  are  provided  in  the  corresponding  input  and  output  data  files. 
These  show  the  use  of  the  code  for  cases  of  (1)  nonuniformly  spaced  points,  and  (2)  several 
heat  elements  with  various  weights. 

These  are  contained  in  the  input2.dat  and  input3.dat  files  with  the  corresponding  output 
data  files  of  output2.dat  and  output3.dat.  The  example  ofinput2.dat  and  output2.dat  is 
annotated  in  tables  5  and  6. 

This  is  for  the  case  of  a  thick  middle  layer  in  the  structure  which  can  be  handled  with 
TXYZ20.  The  example  of  input3.dat  and  output3.dat  is  illustrated  in  tables  7  and  8.  The 
specific  example  shows  how  to  construct  the  problem  for  the  case  of  two  elements  where 
one  is  a  heat  source  and  the  other  is  a  heat  sink.  Nonequispaced  sets  of  x  and  z  values  are 
also  used  in  the  calculation.  These  examples  contain  all  of  the  features  of  the  use  of  the 
TXYZ20  code. 

All  of  the  input  and  output  data  files  are  provided  to  give  the  user  hands-on  experience 
with  the  use  of  the  code  and  to  make  sure  that  the  code  is  properly  working  on  the  user's 
machine. 

Values  of  the  thermal  conductivity  of  materials  of  interest  to  the  semiconductor  community 
may  be  found  in  reference  [5].  The  use  of  the  basic  TXYZ  code  for  thermal  characterization 
is  discussed  there. 

The  code  has  also  been  used  in  the  original  TXYZ  form  for  the  modeling  of  MMIC  devices 
(GaAs)  for  the  determination  of  the  channel  temperature  when  the  devices  are  undergoing 
life  testing  [6]. 

Finally,  the  field  of  thermal  resistance  measurements  has  been  reviewed  [7],  This  reference 
contains  results  which  are  related  to  the  use  of  the  original  TXYZ  code. 
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Table  3.    CPU  times  used  in  TXYZ20  calculation  using  inputl.dat  data. 


This  table  presents  the  CPU  time  needed  to  carry  out  the  calculation 
using  the  data  in  inputl.dat  (annotated  in  table  1)  on  several  machines. 
This  information  is  intended  to  provide  a  benchmark  for  the  running  of 
the  TXZY20  code.    The  machines  represented  were  used  due  to  availability 
and  do  not  represent  a  recommendation.    The  relative  speeds  presented 
cover  a  good  deal  of  the  spectrum  to  be  expected. 

Machine  CPU  Relative  speed 


VAX  457  sec.  1 

SUN  205  sec.  2.2 

SPARC  2  60  sec.  7.6 

CRAY  8.375  sec.  55 


Table  4.    Annotated  input  data  file  (inputla.dat)  for  TXYZ20  code. 

The  data  correspond  to  those  contained  in  the  file  inputla.dat  which  can  be 

used  to  run  the  case  of  a  single  heat  source  of  unit  weight.    This  illustrates 

the  use  of  the  IEDGEX  =  2  generation  of  a  set  of  values  of  x  to  be  used  in  the 

calculation.    The  final  set  of  data  is  equivalent  to  that  of  inputl.dat  for 

IEDGE  =  1.  The  output  is  the  same  as  that  contained  in  outputl.dat  (discussed 
i  n  table  2) . 

INPUT  DATA  DESCRIPTION  OF  INPUT  DATA 

100  100  X  AND  Y  DIMENSIONS  OF  RECTANGULAR  STRUCTURE 

1  .001  THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  TOP  LAYER 

5  .0001  THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  MIDDLE  LAYER 

4  .005  THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  BOTTOM  LAYER 

500  500  UPPER  SUMMATION  LIMITS  FOR  N  AND  M  SUMMATIONS 

2  IEDGEX  (=2,  THEN  READ  IN  INDIVIDUAL  VALUES) 
11  NUMBER  OF  VALUES 

0.0  X(l)      FIRST  VALUE  OF  X 

10.  X(2)      SECOND  VALUE  OF  X 

20.  X(3) 

30.  X(4) 

40.  X(5) 

50.  X(6) 

60.  X(7) 

70.  X(8) 

80.  X(9) 

90.  X(10) 

100.  X(ll)    ELEVENTH  VALUE  OF  X 

1  IEDGEY  (=1,  THEN  READ  IN  THREE  VALUES  ON  NEXT  LINE) 

1    50.    0.0  NUMBER  OF  VALUES,  FIRST  POINT,  INCREMENT 

1  IEDGEZ  (=1,  THEN  READ  IN  THREE  VALUES  ON  NEXT  LINE) 

1    0.0    0.0  NUMBER  OF  VALUES,  FIRST  POINT,  INCREMENT 

1      1.0  NUMBER  OF  SOURCES  AND  POWER  DENSITY 

1    40.0    20.0    40.0    20.0  WEIGHT,  X,  LENGTH  ALONG  X,  Y,  LENGTH  ALONG  Y 
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Table  5.    Annotated  input  data  file  (input2.dat)  for  TXYZ20  code. 


The  data  correspond  to  those  contained  in  the  file  input2.dat  which  is  used 
to  illustrate  the  use  of  the  IEDGEZ  =  2  generation  of  a  set  of  values  of  z 
to  be  used  in  the  calculation.    This  is  for  the  case  of  a  thick  middle  layer 
where  the  original  TXYZ  code  had  difficulty  with  an  artificial  overflow. 
The  output  for  this  case  is  contained  in  output2.dat  and  is  discussed  in 
table  6. 


INPUT  DATA  DESCRIPTION  OF  INPUT  DATA 

1060.00     280.00  X  AND  Y  DIMENSIONS  OF  RECTANGULAR  STRUCTURE 

3.0    .00011  THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  TOP  LAYER 

250    .00011  THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  MIDDLE  LAYER 

0.1    .00011  THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  BOTTOM  LAYER 

50     50  UPPER  SUMMATION  LIMITS  FOR  N  AND  M  SUMMATIONS 

1  IEDGEX  (=1,  THEN  READ  IN  THREE  VALUES  ON  NEXT  LINE) 

10.    0.  NUMBER  OF  VALUES,  FIRST  POINT,  INCREMENT 

1  IEDGEY  (=1,  THEN  READ  IN  THREE  VALUES  ON  NEXT  LINE) 
10.    0.  NUMBER  OF  VALUES,  FIRST  POINT,  INCREMENT 

2  IEDGEZ  (=2,  THEN  READ  IN  INDIVIDUAL  VALUES) 
10  NUMBER  OF  VALUES 

0.  Z(l)      FIRST  VALUE  OF  Z 

1.0  Z(2)      SECOND  VALUE  OF  Z 

2.0  Z(3) 

3.0  Z(4) 

4.0  Z(5) 

4.5  Z(6) 

100.  Z(7) 

200.  Z(8) 

230.  Z(9) 

240.  Z(10)    TENTH  VALUE  OF  Z 

1      1.0  NUMBER  OF  SOURCES  AND  POWER  DENSITY 

1    0.    53.      0.    14.  WEIGHT,  X,  LENGTH  ALONG  X,  Y,  LENGTH  ALONG  Y 
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Table  6.    Annotated  output  data  file  (output2.dat)  for  TXYZ20  program. 

The  data  correspond  to  those  contained  in  the  file  output2.dat  which  is 
generated  from  the  input  data  contained  in  input2.dat  (see  table  5). 
The  first  block  of  lines  are  the  x,y,z,  and  calculated  t(x,y,z)  values. 
Beyond  these  are  the  listing  of  the  input  variables. 


x 

y 

z 

t  ( X   V  . 

0 

0000000E+00 

0 

0000000E+00 

0.0000000E+00 

233445.2 

0 

0000000E+00 

0 

0000000E+00 

-1.000000 

224832.5 

0 

00O00OOE+0O 

0 

0000000E+00 

-2.000000 

216600.1 

0 

0OO00O0E+00 

0 

0000000E+00 

-3.000000 

208748.5 

0 

0000000E+00 

0 

0000000E+00 

-4.000000 

201273.0 

0 

0O0OO00E+O0 

0 

0000000E+00 

-4.500000 

197674.8 

0 

00000O0E+00 

0 

0000000E+00 

-100.0000 

29815.45 

0 

0000000E+00 

0 

0000000E+00 

-200.0000 

7110.651 

0 

0000000E+00 

0 

0000000E+00 

-230.0000 

2989.828 

0 

0000000E+00 

0 

0000000E+00 

-240.0000 

1686.604 

STEADY-STATE  THERMAL  ANALYSIS  CALCULATION  USING  EQS.  (13) -(23)  OF  KOKKAS 

THERMAL  CONDUCTIVITIES  AND  LAYER  THICKNESSES 
LX=  1060.00       LY=  280.00 

Ll=  3.00000  L2=  250.00000  L3=  0.10000 
Kl=  0.00011000  K2=  0.00011000  K3=  0.00011000 

UPPER  SUMMATION  LIMITS       NUP=     50    MUP=  50 


NUMBER  OF  HEAT  S0URCES=  1 
POWER  DENSITY=  1.000000 

COORDINATES,  LENGTHS,  AND  WIDTHS  OF  HEAT  SOURCES 

HEAT  SOURCE       WTSOUR        XSOUR  YSOUR  LXSOUR  LYSOUR 

1  1.00000  0.00000  0.00000         53.00000  14.00000 
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Table  7.    Annotated  input  data  file  (input3.dat)  for  TXYZ20  code. 


The  data  correspond  to  those  contained  in  the  file  input3.dat  which  is  used 
to  illustrate  the  use  of  the  IEDGEX  =  2  and  IEDGEZ  =  2  generation  of  a  set 

of  values  of  x  and  z  to  be  used  in  the  calculation.     In  addition,  two  elements 

are  used.  One  is  a  source  with  a  weight  of  1.0  and  the  other  is  a  sink  with  a 

weight  of  -1.0.    The  output  for  this  case  is  contained  in  output2.dat  and  is 

discussed  in  table  8. 

INPUT  DATA  DESCRIPTION  OF  INPUT  DATA 

100  100  X  AND  Y  DIMENSIONS  OF  RECTANGULAR  STRUCTURE 

1  .001  THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  TOP  LAYER 

5  .0001  THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  MIDDLE  LAYER 

4  .005  THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  BOTTOM  LAYER 

500  500  UPPER  SUMMATION  LIMITS  FOR  N  AND  M  SUMMATIONS 

2  IEDGEX  (=2,  THEN  READ  IN  INDIVIDUAL  VALUES) 
10  NUMBER  OF  VALUES 

40.  X(l)      FIRST  VALUE  OF  X 

45.  X(2)      SECOND  VALUE  OF  X 

46.  X(3) 

47.  X(4) 

48.  X(5) 

52.  X(6) 

53.  X(7) 

54.  X(8) 

55.  X(9) 

60.  X(10)    TENTH  VALUE  OF  X 

1  IEDGEY  (=1,  THEN  READ  IN  THREE  VALUES  ON  NEXT  LINE) 

1  50.    0.0  NUMBER  OF  VALUES,  FIRST  POINT,  INCREMENT 

2  IEDGEZ  (=2,  THEN  READ  IN  INDIVIDUAL  VALUES) 
9  NUMBER  OF  VALUES 

0.  0  Z(l)      FIRST  VALUE  OF  Z 
1.0  Z(2)      SECOND  VALUE  OF  Z 
2.0  Z(3) 

3.0  Z(4) 

4.0  Z(5) 

5.0  Z(6) 

6.0  Z(7) 

9.0  Z(8) 

10.0  Z(9)      NINTH  VALUE  OF  Z 

2     1.0  NUMBER  OF  SOURCES  AND  POWER  DENSITY 

1.  35.0  10.0  45.0  10.0  WEIGHT,  X,  LENGTH  ALONG  X,  Y,  LENGTH  ALONG  Y 
-1.    55.0  10.0  45.0    10.0    WEIGHT,  X,  LENGTH  ALONG  X,  Y,  LENGTH  ALONG  Y 
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Table  8.    Annotated  output  data  file  (output3.dat)  for  TXYZ20  program. 

The  data  correspond  to  those  contained  in  the  file  output3.dat  which  i 
generated  from  the  input  data  contained  in  input3.dat  (see  table  7). 
The  first  block  of  lines  are  the  x,y,z,  and  calculated  t(x,y,z)  values 
Beyond  these  are  the  listing  of  the  input  variables.    This  example  has 
been  constructed  to  show  not  only  variable  x  and  z  but  also  the  use  of 
negative  weights  on  the  elements  in  the  construction  of  heat  sinks.  I 
this  particular  case,  the  two  elements  are  of  equal  size  are  are  sym- 
metrical ly  placed  with  respect  to  the  middle  of  the  structure.  The 
calculated  t(x,y,z)  values  directly  reflect  the  opposite  but  equal 
weights  of  the  two  elements. 


x 

y 

z 

t(x,y,z) 

40 

00000 

50 

00000 

0. 

OOOOOOOE+00 

13895.94 

45 

00000 

50 

00000 

0. 

OOOOOOOE+00 

8401. 

479 

46 

00000 

50 

00000 

0. 

OOOOOOOE+00 

6187. 

981 

47 

00000 

50 

00000 

0. 

OOOOOOOE+00 

4418. 

474 

48 

00000 

50 

00000 

0. 

OOOOOOOE+00 

2844. 

606 

52 

00000 

50 

00000 

0. 

OOOOOOOE+00 

-2844. 

609 

53 

00000 

50 

00000 

0. 

OOOOOOOE+00 

-4418. 

476 

54 

00000 

50 

00000 

0. 

OOOOOOOE+00 

-6187. 

989 

55 

00000 

50 

00000 

0. 

OOOOOOOE+00 

-8401 . 

491 

60 

00000 

50 

00000 

0. 

OOOOOOOE+00 

-13895.94 

40 

00000 

50 

00000 

-1 

.000000 

13214.03 

45 

00000 

50 

00000 

-1 

.000000 

8048. 

078 

46 

00000 

50 

00000 

-1 

.000000 

6112. 

896 

47 

00000 

50 

00000 

-1 

.000000 

4373. 

840 

48 

00000 

50 

00000 

-1 

.000000 

2817. 

945 

52 

00000 

50 

00000 

-1 

.000000 

-2817. 

947 

53 

00000 

50 

00000 

-1 

.000000 

-4373. 

844 

54 

00000 

50 

00000 

-1 

.000000 

-6112. 

901 

55 

00000 

50 

00000 

-1 

.000000 

-8048 . 

083 

60 

00000 

50 

00000 

-1 

.000000 

-13214.03 

40 

00000 

50 

00000 

-2 

.000000 

9877. 

484 

45 

00000 

50 

00000 

-2 

.000000 

6111. 

877 

46 

00000 

50 

00000 

-2 

.000000 

4777. 

800 

47 

00000 

50 

00000 

-2 

.000000 

3481. 

801 

48 

00000 

50 

00000 

-2 

.000000 

2265. 

086 

52 

00000 

50 

00000 

-2 

.000000 

-2265. 

088 

53 

00000 

50 

00000 

-2 

.000000 

-3481. 

803 

54 

00000 

50 

00000 

-2 

.000000 

-4777. 

805 

55 

00000 

50 

00000 

-2 

.000000 

-6111. 

882 

60 

00000 

50 

00000 

-2 

.000000 

-9877. 

482 

40 

00000 

50 

00000 

-3 

.000000 

7033. 

427 

45 

00000 

50 

00000 

-3 

.000000 

4407. 

094 

46 

00000 

50 

00000 

-3 

.000000 

3498. 

253 

47 

00000 

50 

00000 

-3 

.000000 

2583. 

641 

48 

.00000 

50 

00000 

-3 

.000000 

1695. 

754 

52 

.00000 

50 

.00000 

-3 

.000000 

-1695. 

755 

53 

.00000 

50 

.00000 

-3 

.000000 

-2583. 

643 

54 

.00000 

50 

.00000 

-3 

.000000 

-3498 

257 

55 

.00000 

50 

.00000 

-3 

.000000 

-4407 . 

097 

60 

.00000 

50 

.00000 

-3 

.000000 

-7033 

424 

15 


X 

40.00000 
45.00000 
46.00000 
47.00000 
48.00000 
52.00000 
53.00000 
54.00000 
55.00000 
60.00000 
40.00000 
45.00000 
46.00000 
47.00000 
48.00000 
52.00000 
53.00000 
54.00000 
55.00000 
60.00000 
40.00000 
45.00000 
46.00000 
47.00000 
48.00000 
52.00000 
53.00000 
54.00000 
55.00000 
60.00000 
40.00000 
45.00000 
46.00000 
47.00000 
48.00000 
52.00000 
53.00000 
54.00000 
55.00000 
60.00000 
40.00000 
45.00000 
46.00000 
47.00000 
48.00000 
52.00000 
53.00000 
54.00000 
55.00000 
60.00000 


y 

50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 
50.00000 


-4.000000 
-4.000000 
-4.000000 
-4.000000 
-4.000000 
-4.000000 
-4.000000 
-4.000000 
-4.000000 
-4.000000 
-5.000000 
-5.000000 
-5.000000 
-5.000000 
-5.000000 
-5.000000 
-5.000000 
-5.000000 
-5.000000 
-5.000000 
-6.000000 
-6.000000 
-6.000000 
-6.000000 
-6.000000 
-6.000000 
-6.000000 
-6.000000 
-6.000000 
-6.000000 
-9.000000 
-9.000000 
-9.000000 
-9.000000 
-9.000000 
-9.000000 
-9.000000 
-9.000000 
-9.000000 
-9.000000 
-10.00000 
-10.00000 
-10.00000 
-10.00000 
-10.00000 
-10.00000 
-10.00000 
-10.00000 
-10.00000 
-10.00000 


t(x,y,z) 

4544.860 

2873.673 

2300.854 

1713.657 

1131.903 
-1131.905 
-1713.658 
-2300.855 
-2873.675 
-4544.859 

2285.325 

1453.400 

1169.267 

875.0942 

580.2814 
-580.2822 
-875.0949 
-1169.268 
-1453.401 
-2285.325 

139.8420 

90.59629 

73.84338 

56.01737 

37.57924 
-37.57928 
-56.01740 
-73.84344 
-90.59635 
-139.8419 

31.62482 

20.66124 

16.93877 

12.92600 

8.715524 
-8.715535 
-12.92601 
-16.93879 
-20.66125 
-31.62481 
0.0000000E+00 
0.0000000E+00 
0.0000000E+00 
0.0000000E+00 
0.0000000E+00 
0.0000000E+00 
0.0000000E+00 
0.0000000E+00 
0.0000000E+00 
0.0000000E+00 
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STEADY-STATE  THERMAL  ANALYSIS  CALCULATION  USING  EQS.  (13) -(23)  OF  KOKKAS 


THERMAL  CONDUCTIVITIES  AND  LAYER  THICKNESSES 
LX=    100.00       LY=  100.00 

Ll=  1.00000  L2=  5.00000  L3=  4.00000 
Kl=  0.00100000  K2=  0.00010000  K3=  0.00500000 

UPPER  SUMMATION  LIMITS       NUP=    500    MUP=  500 


NUMBER  OF  HEAT  SOURCES=  2 
POWER  DENSITY=  1.000000 

COORDINATES,  LENGTHS,  AND  WIDTHS  OF  HEAT  SOURCES 
HEAT  SOURCE       WTSOUR        XSOUR  YSOUR 


1 

2 


1.00000 
-1.00000 


35.00000 
55.00000 


45.00000 
45.00000 


LXSOUR 

10.00000 
10.00000 


LYSOUR 

10.00000 
10.00000 
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APPENDIX 


C*******************************************^ 

C  TXYZ  VERSION  2.0 

C 

C  THIS  PROGRAM  CALCULATES  THE  STEADY-STATE  TEMPERATURE  DISTRIBUTION, 

C  T(X,Y,Z),  DUE  TO  DC  POWER  INPUTS/OUTPUTS.    THE  TEMPERATURE  IS 

C  RELATIVE  TO  THE  AMBIENT.    THE  SPECIFIC  EQUATIONS  USED  ARE  GIVEN  IN 

C  EQUATIONS  (13) -(23),  WITH  S=0  (ZERO  FREQUENCY,  STEADY-STATE  CONDITION), 

C  IN  THE  PAPER  BY  KOKKAS. 

C 

C  REFERENCES:  THE  ORIGINAL  MATHEMATICAL  ANALYSIS  OF  THE  THREE  LAYER 

C  STRUCTURE  WAS  PERFORMED  IN    THE  PAPER  "THERMAL  ANALYSIS 

C  OF  MULTIPLE-LAYERED  STRUCTURES"  BY  ACHILLES  G.  KOKKAS, 

C  IEEE  TRANS.  ELEC.  DEV.  VOL.  ED-21,  NO.  11,  674-681  (1974)). 

C  THE  ORIGINAL  FORTRAN  IMPLEMENTATION  OF  THE  STEADY  STATE 

C  IN  THE  FORM  OF  THE  ORIGINAL  TXYZ  CODE  IS  PRESENTED  IN  THE 

C  REPORT  "SEMICONDUCTOR  MEASUREMENT  TECHNOLOGY:  TXYZ:  A 

C  PROGRAM  FOR  SEMICONDUCTOR  IC  THERMAL  ANALYSIS"    BY  JOHN 

C  ALBERS,  NBS  SPEC.  PUBL.  400-76  (APRIL  1984) 

C 

C  THE  CALCULATION  IS  PERFORMED  FOR  A  RECTANGULAR  THREE  LAYER  STRUCTURE 

C  WITH  AN  ARBITRARY  NUMBER  OF  RECTANGULAR  HEAT  SOURCES/SINKS  ON  THE  TOP 

C  SURFACE.    THE  TEMPERATURE  MAY  BE  CALCULATED  AT  ANY  POINT  (X,Y,Z)  INSIDE 

C  OR  ON  THE  SURFACES  OF  THE  STRUCTURE.    THE  CALCULATION  FOLLOWS  FROM  THE 

C  INPUT  OF  THE  THICKNESSES  AND  THERMAL  CONDUCTIVITIES  OF  THE  THREE  LAYERS. 

C  IT  IS  IMPORTANT  TO  EMPHASIZE  THAT  THE  CALCULATION  IS  GENERAL  FOR  THE 

C  THREE  LAYER  STRUCTURE  AND  THE  APPLICATION  TO  SEMICONDUCTOR  STRUCTURES 

C  IS  A  SPECIAL  CASE. 
C 

C***************************************^ 
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1 

2 

3 
4 
5 
6 

7 
8 
9 

10 
11 
12 

13 
14 

15 
16 

17 
27 
22 
31 
51 
52 
53 
54 
55 
56 
57 
88 


DIMENSION  X(100),  Y(100),  Z(100),  COSYT(500) 

DIMENSION  ARUZER(500,500),  ARFUNZ(500,500) 

DIMENSION  WTS0UR(5O),  XS0UR(5O),  YS0UR(5O), 

REAL  LXSOUR(50),  LYSOUR(50),  Kl,  K2,  K3,  LX,  LY,  LI,  L2,  L3 

COMMON    Kl,  K2,  K3,  LX,  LY,  LI,  L2,  L3 

COMMON  NSOUR,WTSOUR,XSOUR,YSOUR,  LXSOUR,  LYSOUR 

PI=3. 14159265 


FORMAT 
FORMAT 
1  (13)- 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


MUP=',I5/) 


1H1) 

IX, 'STEADY-STATE  THERMAL  ANALYSIS  CALCULATION  USING  EQS. 
23)  OF  KQKKAS'/) 

IX, 'THERMAL  CONDUCTIVITIES  AND  LAYER  THICKNESSES') 
IX, 'Kl=  ',F10.8,'  K2=  ',F10.8,'  K3=  ',F10.8) 
1X,'L1=  ',F10.5,'  L2=  ',F10.5,'  L3=  ',F10.5) 
/IX, 'UPPER  SUMMATION  LIMITS  ',2X,'  NUP=',I5, 


//IX, 'NUMBER  OF  HEAT  SOURCES=' ,15) 

/IX, 'COORDINATES,  LENGTHS,  AND  WIDTHS  OF  HEAT  SOURCES'/) 
IX, 'HEAT  SOURCE  ' ,3X, 'WTSOUR' ,5X, 'XSOUR' ,9X, 'YSOUR' ,8X, 


FORMAT 

FORMAT 

FORMAT 
1  ' LXSOUR ',8X, 'LYSOUR'/) 

FORMAT 

FORMAT 

FORMAT 
1  ,F11 

FORMAT 

FORMAT 
1  ,F11 

FORMAT 

FORMAT 
1  ,F11 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


4X,I3,5X,F10.5,3X,F10.5,3X,F10.5,3X,F10.5,3X,F10.5) 
IX, 'POWER  DENSITY=',F11.6) 

/IX, 'CALCULATING  ',13,'  X  POINTS  WITH  A  FIRST  POINT  OF 
6,'  AND  A  STEP  SIZE  OF  ',F11.6) 
/IX, 'THE  CONSTANT  X  COORDINATE  IS  \F11.6) 
/IX, 'CALCULATING  ',13,'  Y  POINTS  WITH  A  FIRST  POINT  OF 
6,'  AND  A  STEP  SIZE  OF  ',F11.6) 
/IX, 'THE  CONSTANT  Y  COORDINATE  IS  ',F11.6) 
/IX, 'CALCULATING  ',13,'  Z  POINTS  WITH  A  FIRST  POINT  OF 
6,'  AND  A  STEP  SIZE  OF  ',F11.6) 
/IX, 'THE  CONSTANT  Z  COORDINATE  IS  ',F11.6) 
1X,'LX=  ',F7.2,3X,'  LY=  ',F7.2) 
1X,F12.4,2X,F12.4,2X,F12.4,2X,F12.4) 
IX, 617) 
IX, 14, 3X, 14) 
1X,F10.5,3X,F10.5) 
IX, II) 

1X,I4,3X,F10.5,3X,F10.5) 
1X,F10.5) 
1X,I2,3X,F10.5) 

1X,F10.5,3X,F10.5,3X,F10.5,3X,F10.5) 


FORMAT (IX, 'YOUR  UPPER  LIMIT  OF  SUMMATION  IS  TOO  LARGE.  TRY  AGAIN') 
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c*************************************^ 
c 

C   INPUT  SECTION  

C        THE  FOLLOWING  VARIABLES  ARE  READ  BY  THE  PROGRAM  FROM  FOROIO 
C 

C  LX  X  DIMENSION  OF  3  LAYER  STRUCTURE 

C  LY  Y  DIMENSION  OF  3  LAYER  STRUCTURE 

C  LI  THICKNESS  OF  TOP  LAYER 

C  Kl  THERMAL  CONDUCTIVITY  OF  TOP  LAYER 

C  L2  THICKNESS  OF  MIDDLE  LAYER 

C  K2  THERMAL  CONDUCTIVITY  OF  MIDDLE  LAYER 

C  L3  THICKNESS  OF  BOTTOM  LAYER 

C  K3  THERMAL  CONDUCTIVITY  OF  BOTTOM  LAYER 

C  NUP  UPPER  LIMIT  OF  N  SUM,  X  DIRECTION 

C  MUP  UPPER  LIMIT  OF  M  SUM,  Y  DIRECTION 

C 

C  IEDGEX  INDEX  FOR  GENERATING  THE  VALUES  OF  X  TO  BE  USED 

C  =1  TO  READ  DATA  FOR  FIXED  INCREMENT  X  VALUES 

C  =2  TO  READ  IN  ARRAY  OF  X  VALUES  OF  NONFIXED  INCREMENT) 

C  IF  IEDGEX=1  THEN  READ  THE  THREE  VARIABLES  (ON  SAME  LINE) 

C  ILX    THE  NUMBER  OF  X  VALUES  TO  BE  USED 

C  XI     THE  VALUE  OF  THE  FIRST  POINT  IN  X 

C  STEPX  (THE  INCREMENT  IN  X) 

C  IF  IEDGEX=2  THEN  READ  THE  VARIABLE  AND  ARRAY  (ONE  PER  LINE) 

C  ILX    THE  NUMBER  OF  X  VALUES  TO  BE  USED 

C  X(I)  THE  ARRAY  OF  X  VALUES  (1=1, ILX) 
C 

C  IEDGEY  INDEX  FOR  GENERATING  THE  VALUES  OF  Y  TO  BE  USED 

C  =1  TO  READ  DATA  FOR  FIXED  INCREMENT  Y  VALUES 

C  =2  TO  READ  IN  ARRAY  OF  X  VALUES  OF  NONFIXED  INCREMENT) 

C  IF  IEDGEY=1  THEN  READ  THE  THREE  VARIABLES  (ON  SAME  LINE) 

C  ILY    THE  NUMBER  OF  Y  VALUES  TO  BE  USED 

C  Yl     THE  VALUE  OF  THE  FIRST  POINT  IN  Y 

C  STEPY  (THE  INCREMENT  IN  Y) 

C  IF  IEDGEY=2  THEN  READ  THE  VARIABLE  AND  ARRAY  (ONE  PER  LINE) 

C  ILY    THE  NUMBER  OF  Y  VALUES  TO  BE  USED 

C  Y(I)  THE  ARRAY  OF  Y  VALUES  (1=1, ILY) 
C 

C  IEDGEZ  INDEX  FOR  GENERATING  THE  VALUES  OF  Z  TO  BE  USED 

C  =1  TO  READ  DATA  FOR  FIXED  INCREMENT  Z  VALUES 

C  =2  TO  READ  IN  ARRAY  OF  Z  VALUES  OF  NONFIXED  INCREMENT) 

C  IF  IEDGEZ=1  THEN  READ  THE  THREE  VARIABLES  (ON  SAME  LINE) 

C  ILZ    THE  NUMBER  OF  Z  VALUES  TO  BE  USED 

C  Zl     THE  VALUE  OF  THE  FIRST  POINT  IN  Z 

C  STEPZ  (THE  INCREMENT  IN  Z) 

C  IF  IEDGEZ=2  THEN  READ  THE  VARIABLE  AND  ARRAY  (ONE  PER  LINE) 

C  ILZ    THE  NUMBER  OF  Z  VALUES  TO  BE  USED 

C  Z(I)  THE  ARRAY  OF  Z  VALUES  (1=1, ILZ) 

C 

C  NOTE:  ENTER  THE  Z-RELATED  VARIABLES  (Zl,  STEPZ  OR  THE  Z(I)  ARRAY) 

C  AS  ZERO  OR  POSITIVE  QUANTITIES.    THE  PROGRAM  CONVERTS  THE  FINAL 

C  Z(I)  ARRAY  TO  ZERO  OR  NEGATIVE  QUANTITIES  AS  THE  CALCULATION 

C  TAKES  THE  Z  VARIABLE  TO  BE  ZERO  OR  NEGATIVE. 

C 
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C  NSOUR  NUMBER  OF  HEAT  SOURCES  (UP  TO  50) 

C  po  POWER  DENSITY 

C 

C  THE  NEXT  NSOUR  LINES  READ  THE  FOLLOWING  INFORMATION  FOR  THE 

C  HEAT  SOURCES  (WITH  ALL  THE  INFORMATION  FOR  EACH  OF  THE  ELEMENTS 

C  ON  A  SINGLE  LINE) 

C 

C  WTSOUR (I) -WEIGHTING  FACTOR  OF  I-TH  SOURCE 

C  (POSITIVE  FOR  SOURCE,  NEGATIVE  FOR  SINK) 

C  XSOUR(I)— X  COORDINATE  OF  ORIGIN  OF  I-TH  SOURCE 

C  LXSOUR (I) -LENGTH  ALONG  X  DIRECTION  OF  I-TH  SOURCE 

C  YSOUR(I) — Y  COORDINATE  OF  ORIGIN  OF  I-TH  SOURCE 

C  LYSOUR (I) -LENGTH  ALONG  Y  DIRECTION  OF  I-TH  SOURCE 

C 

C****************************************^ 

C         READ  LX  AND  LY  (THE  X  AND  Y  DIMENSIONS  OF  THE  RECTANGULAR  STRUCTURE) 
READ (10,*)  LX,  LY 

C        READ  LI  AND  Kl  (THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  THE  TOP  LAYER) 
READ(10,*)  LI,  Kl 

C        READ  L2  AND  K2  (THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  THE  MIDDLE  LAYER) 
READ (10,*)  L2,  K2 

C        READ  L3  AND  K3  (THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  THE  BOTTOM  LAYER) 
READ (10,*)  L3,  K3 

C        READ  NUP  AND  MUP  (UPPER  LIMIT  OF  THE  SUMMATION  OVER  THE  INDEX  N  (X-DIR) 
C  UPPER  LIMIT  OF  THE  SUMMATION  OVER  THE  INDEX  M  (Y-DIR)) 

READ (10,*) NUP, MUP 

C         NUP  AND  MUP  MUST  BE  LESS  THAN  OR  EQUAL  TO  THE  DIMENSIONALITY  OF 
IF  (NUP. GT. 500. OR. MUP. GT. 500)  GO  TO  3999 
READ(10,*)IEDGEX 
GOTO  (110,115)IEDGEX 


110  READ  (10,*)ILX,X1,STEPX 
DO  111  1=1, ILX 
X(I)=X1+(I-1)*STEPX 

111  CONTINUE 
GOTO  119 

115  READ  (10,*) ILX 
DO  116  1=1, ILX 
READ(10,*)X(I) 

116  CONTINUE 

119  READ(10,*)IEDGEY 
GOTO  (120,125)  IEDGEY 

120  READ  (10,*)ILY,Y1,STEPY 
DO  121  1=1, ILY 
Y(I)=Y1+(I-1)*STEPY 

121  CONTINUE 
GOTO  129 

125  READ  (10,*) ILY 
DO  126  1=1, ILY 
READ(10,*)Y(I) 

126  CONTINUE 
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129  READ(10,*)IEDGEZ 
GOTO  (130,135)  IEDGEZ 

130  READ  (10,*)ILZ,Z1,STEPZ 
Z1=-1.0*Z1 
STEPZ=-1.0*STEPZ 

DO  131  1=1, ILZ 
Z(I)=Z1+(I-1)*STEPZ 

131  CONTINUE 
GOTO  139 

135  READ  (10,*) ILZ 
DO  136  1=1, ILZ 
READ(10,*)Z(I) 
Z(I)  =  -1.0*Z(I) 

136  CONTINUE 

C        READ  THE  NUMBER  OF  HEAT  SOURCES  AND  THE  POWER  DENSITY 

C        NOTE-POWER  DENSITY  IS  MULTIPLICATIVE  FACTOR  USUALLY  SET  EQUAL  TO  UNITY 

139    READ (10,*) NSOUR, PO 
C        PO  IS  THE  POWER  DENSITY,  ASSUMED  UNIFORM  FOR  ALL  HEATERS 
C        NSOUR  IS  THE  TOTAL  NUMBER  OF  HEATING  ELEMENTS  ON  THE  SURFACE  OF  THE 
C        THE  TOP  LAYER  (UP  TO  50) 

C        THE  NEXT  LOOP  READS  IN  THE  COORDINATES  OF  THE  ORIGIN  OF  THE 
C        HEATING  ELEMENTS  ALONG  WITH  THEIR  LENGTHS  AND  WIDTHS 
C        THE  WEIGHTING  FACTOR  IS  ALSO  ENTERED  (THIS  IS  REAL,  NONINTEGER) 
DO  100  1=1, NSOUR 

READ(10,*)WTSOUR(I),XSOUR(I),  LXSOUR(I) ,  YSOUR(I)  ,  LYSOUR(I) 
C        WTSOUR(I)  IS  THE  WEIGHTING  FACTOR  FOR  THE  I-TH  HEATER  ELEMENT 
C        XSOUR(I)  IS  THE  X  COORDINATE  OF  THE  ORIGIN  OF  I-TH  HEATER  ELEMENT 
C        LXSOUR(I)  IS  THE  LENGTH  OF  THE  I-TH  HEATER  ALONG  THE  X  DIRECTION 
C        YSOUR(I)  IS  THE  Y  COORDINATE  OF  THE  ORIGIN  OF  I-TH  HEATER  ELEMENT 
C        LYSOUR(I)  IS  THE  LENGTH  OF  THE  I-TH  HEATER  ALONG  THE  Y  DIRECTION 
100  CONTINUE 

P04LK  =  4.0  *  PO  /  (  LX*  LY*  Kl) 

PILX  =  PI  /  LX 

PILY  =  PI  /  LY 

C*************************************************************************** 

c 

C         END  OF  DATA  INPUT  SECTION 

C        BEGIN  CALCULATION  OF  T(X,Y,Z) 

C        THE  SUBROUTINES  USED  IN  THE  CALCULATION  ARE: 

C         1)  UZERO(N,M)  -  CALCULATES  THE  FOURIER  COSINE  TRANSFORM  OF  THE 

C  FUNCTION,  U(X,Y),  THE  POWER  DENSITY  FUNCTION  FOR  ALL  OF  THE 

C  HEAT  SOURCES. 

C        2)  FUNZ(N,M,Z)  -  CALCULATES  THE  Z-DEPENDENT  PORTION  OF  THE  SUM 
C  REMEMBERING  THAT  THIS  IS  A  FUNCTION  OF  THE  SUMMATION 

C  INDICES  (N,M). 

p*************************************************************************** 
C         CALCULATE  THE  FOURIER  COMPONENTS  OF  THE  HEAT  SOURCES,  U(N,M) 
C********* ************************************************** 

DO  300  MM=1,MUP 

M  =  MM  -  1 

DO  250  NN=1,NUP 
N  =  NN  -  1 

ARUZER (NN , MM) =UZERO (N , M) 
250  CONTINUE 
300  CONTINUE 
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c*****************************************^ 
c 

C        END  OF  U(N,M)  CALCULATION  AND  BEGINNING  OF  MAJOR  LOOP  FOR  Z 
C 

C*****************************************^ 
DO  3000  IZ=1,ILZ 

c***************************************^ 

C        CALCULATE  THE  Z  DEPENDENT  POTION,  I.E.,  FUNZ(N,M,Z) 
C*******^*********************************^ 

DO  400  MM=1,MUP 

M  =  MM  -  1 

DO  350  NN=1,NUP 
N  =  NN  -  1 

ARFUNZ (NN , MM) =FUNZ  (N , M , Z  (IZ) ) * ARUZER (NN , MM) 
350  CONTINUE 
400  CONTINUE 

DO  3000  IY=1,ILY 
DO  700  MM=1,MUP 
M  =  MM  -  1 

COSYT (MM) =COS (FLOAT (M) *Y (IY) *PILY) 
700  CONTINUE 

DO  3000  IX=1,ILX 
SUM=0.0 

DO  1900  MM=1,MUP 

M  =  MM  -  1 

DO  1700  NN=1,NUP 

N  =  NN  -  1 

NDN=0 

NDM=0 

IF  (N.EQ.O)  NDN=1 
IF  (M.EQ.O)  NDM=1 

TOP  =  ARFUNZ (NN, MM)  *  COS (FLOAT (N) *X (IX) *PILX)  *  COSYT (MM) 
B0TT0M=(NDN+1)*(NDM+1) 
TSUM=T0P/B0TT0M 
SUM=SUM+TSUM 
1700  CONTINUE 
1900  CONTINUE 

TEMP  =  P04LK  *  SUM 
WRITE  (12 , *) X (IX) , Y (IY) , Z (IZ) , TEMP 
C  WRITE  (12,22)X(IX),Y(IY),Z(IZ),TEMP 

3000  CONTINUE 
C  WRITE (12,1) 

WRITE (12, 2) 
WRITE (12, 3) 
WRITE (12, 27)  LX,  LY 
WRITE (12, 5)  LI,  L2,  L3 
WRITE(12,4)  Kl,  K2,  K3 
WRITE(12,6)NUP,MUP 
WRITE(12,7)NS0UR 
WRITE(12,11)  PO 
WRITE (12, 8) 
WRITE (12, 9) 
DO  3888  I=1,NS0UR 

WRITE(12,10)I,WTSOUR(I),XSOUR(I),YSOUR(I),  LXSOUR(I),  LYSOUR(I) 
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3888  CONTINUE 
GO  TO  4000 

3999  WRITE (6, 88) 

4000  STOP 
END 

C****************************************************************************** 
C        END  OF  THE  MAIN  PROGRAM 
C        FUNCTION  UZERO  LISTING 

C****************************************************************************** 
C        DESCRIPTION  OF  THE  FUNCTION  UZERO (N,M) 

C        THIS  FUNCTION  CALCULATES  THE  DOUBLE  FOURIER  COSINE  TRANSFORM 
C        OF  THE  POWER  DENSITY  FUNCTION,  U(X,Y).    THIS  IS  THE  TRANSFORM 
C        FOR  ALL  OF  THE  HEAT  SOURCES.    THE  ASSUMPTION  IS  MADE  THAT  THE 
C        POWER  DENSITY  IS  UNIFORM  AND  EQUAL  TO  UNITY  OVER  THE  SURFACE 
C         OF  THE  HEATING  ELEMENTS.    THAT  IS, 
C  U(X,Y)=1    (XSOUR(I)<=X<=XSOUR(I)+LXSOUR(I)  AND 

C  YSOUR(I)<=Y<=YXOUR(I)+LYSOUR(I)  ). 

C  U(X,Y)=0  OTHERWISE. 

C        UNDER  THESE  CONDITIONS,  IT  IS  POSSIBLE  TO  ANALYTICALLY  EVALUATE 
C        THE  DOUBLE  INTEGRAL  FOR  EACH  HEATING  ELEMENT.    AS  THE  HEATING 
C        ELEMENTS  ARE  ASSUMED  TO  BE  INDEPENDENT,  THE  CONTRIBUTION  FROM 
C        EACH  ELEMENT  MAY  BE  ADDED  TO  OBTAIN  THE  U(N,M)  FOR  ALL. 

C 

C  NONUNIFORM  POWER  DENSITIES  MAY  BE  TAKEN  CARE  OF  BY  USING  THE 
C  WEIGHTING  FACTOR  FOR  EACH  ELEMENT  (READ  IN  THE  MAIN  PROGRAM) 
C****************************************************************************** 

FUNCTION  UZERO (N,M) 

DIMENSION  WTS0UR(5O),  XSOUR(50),  YSOUR(SO), 

REAL  LXS0UR(5O),  LYS0UR(5O),  Kl,  K2,  K3,  LX,  LY,  LI,  L2,  L3 

COMMON    Kl,  K2,  K3,  LX,  LY,  LI,  L2,  L3 

COMMON  NSOUR,WTSOUR,XSOUR,YSOUR,  LXSOUR,  LYSOUR 

PI=3. 14159265 

UZERO=0.0 

DO  500  I=1,NS0UR 

IF(N.EQ.O)  GO  TO  100 

TERMX  =  SIN(FLOAT(N)*PI*(XSOUR(I)  +    LXSOUR (I))/  LX) 
1  -  SIN(FLOAT(N)*PI*XSOUR(I)/  LX) 

TERMX=TERMX*  LX/ (FLOAT (N) *PI) 
GO  TO  150 
100        TERMX=  LXSOUR (I) 
150        IF(M.EQ.O)  GO  TO  200 

TERMY  =  SIN(FLOAT(M)*PI*(YSOUR(I)  +    LYSOUR  (I))/  LY) 
1  -  SIN(FLOAT(M)*PI*YSOUR(I)/  LY) 

TERMY=TERMY*  LY/  (FLOAT (M) *PI) 
GO  TO  250 
200        TERMY=  LYSOUR (I) 
250  TERMI=TERMX*TERMY 

UZERO=UZERO+TERMI*WTSOUR (I) 
500  CONTINUE 
RETURN 
END 

C* ***************************************************************************** 
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C        FUNCTION  FUNZ  LISTING 
C 

C        DESCRIPTION  OF  THE  FUNCTION  FUNZ(N,M,Z) 

C        THIS  FUNCTION  IS  USED  TO  CALCULATE  THE  Z  DEPENDENT  PART  OF  THE 
C        FUNCTION  USING  THE  S=0  VERSIONS  OF  EQUATIONS  (15) -(17)  OF  KOKKAS. 
C        THESE  ARE  USED  IN  CONJUNCTION  WITH  EQUATIONS  (18) -(22)  ALSO  FOR 
C        S=0  (STEADY-STATE  CONDITION) .    THE  SPECIFIC  FORM  OF  FUNZ  IS 
C        DETERMINED  BY  THE  VALUE  OF  Z,  I.E.,  IF  Z  FALLS  IN  THE  TOP,  MIDDLE 
C        OR  BOTTOM  LAYER  OF  THE  STRUCTURE.    IN  ADDITION,  SPECIAL  CARE  IS 
C        TAKEN  AS  TO  EVALUATE  FUNZ  FOR  THE  CASE  WHERE  GAMMA=0  AS  THESE  ARE 
C        SIMPLE,  BUT  THE  COMPUTER  DOES  NOT  KNOW  HOW  TO  EVALUATE  LIMITS. 
C         IN  ADDITION,  THE  FORM  OF  THE  THREE  SOLUTIONS  HAS  BEEN  CHANGED  TO 
C        GET  RID  OF  THE  ARTIFICIAL  OVERFLOW  PROBLEMS  COMING  FROM  THE 
C        HYPERBOLIC  FUNCTIONS,  COSH  AND  SINH,  FOR  THE  CASES  WHERE  THE 
C        ARGUMENTS  BECOME  LARGE. 
C 

C**************************************** 
FUNCTION  FUNZ(N,M,Z) 

DIMENSION  WTSOUR(50),  XSOUR(50),  YSOUR(50), 

REAL  LXSOUR(50),  LYSOUR(50),  Kl,  K2,  K3,  LX,  LY,  LI,  L2,  L3 

COMMON    Kl,  K2,  K3,  LX,  LY,  LI,  L2,  L3 

COMMON  NSOUR,WTSOUR,XSOUR,YSOUR,  LXSOUR,  LYSOUR 

PI=3. 14159265 

GAMMA=SQRT ( (FLOAT (N) *PI/  LX)**2  +  (FLOAT (M) *PI/  LY)**2) 

VS=GAMMA*  LI 

VC=GAMMA*  L2 

VI=GAMMA*  L3 

VT=GAMMA* (  Ll+Z) 

VM=GAMMA* (  L1+  L2+Z) 

VB=GAMMA* (  L1+  L2+  L3+Z) 

BOTl=TANH (VS) *TANH (VI) 

BOTl=BOTl+(  K3/  K2) *TANH (VS) *TANH (VC) 

BOT2=(  K2/  Kl)*TANH(VI)*TANH(VC)+(  K3/  Kl) 

GFUNC=1.0/ (BOT1+BOT2) 

AZ=ABS(Z) 

IF  (AZ.GT.  LI)  GO  TO  500 
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c* ********************************************************** 
C  TOP  LAYER  CALCULATION 

C        THIS  PORTION  IS  THE  TOP  LAYER  CALCULATION  WHICH  IS  DEFAULTED 
C        TO  IF  Z  FALLS  INTO  THE  TOP  LAYER 

C* ***************************************************************************** 
IF  (GAMMA. EQ. 0.0)  GO  TO  100 
TERMSl=TANH(VI)+(  K3/  K2)*TANH(VC) 

TERMS2=(  K2/  Kl)*TANH(VT)*(TANH(VI)*TANH(VC)  +  (  K3/  K2)) 

TERMS=TERMS1+TERMS2 

IF  (Z.EQ.O.O)  GO  TO  90 

IF  (VS.GT.5.0.AND.VT.GT.5.0)  GO  TO  80 

IF  (VS.LT.5.0)  GO  TO  10 

C1=2.0*EXP(-VS) 

GO  TO  20 
10  C1=1.0/COSH(VS) 
20  CONTINUE 

IF  (VT.LT.5.0)  GO  TO  30 

C2=0.5*EXP (VT) 

GO  TO  40 
30  C2=C0SH(VT) 
40  CONTINUE 

FUNZ=GFUNC*TERMS*C1*C2/GAMMA 

RETURN 

80     FUNZ=GFUNC*TERMS*EXP (GAMMA*Z) /GAMMA 
RETURN 

90  FUNZ=GFUNC*TERMS/GAMMA 
RETURN 

100    FUNZ=(  Ll+Z)+(  Kl/  K2)*  L2+(  Kl/  K3)*  L3 

RETURN 
500    TOTAL=  L1+  L2 

IF  (AZ.GT. TOTAL)  GO  TO  1500 
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c******************************************^ 

C  MIDDLE  LAYER  CALCULATION 

C        THIS  IS  THE  MIDDLE  LAYER  CALCULATION  WHICH  IS  DEFAULTED  TO  IF 
C        Z  FALLS  INTO  THIS  DOMAIN  OF  DEPTHS 
C************************************^ 

IF  (GAMMA. EQ. 0.0)  GO  TO  1000 

TERMC=TANH(VI)+(  K3/  K2)*TANH(VM) 

IF  (VS.GT.5.0.AND.VC.GT.5.0.AND.VM.GT.5.0)  GO  TO  800 

IF  (VS.LT.5.0.AND.VC.GT.5.0.AND.VM.GT.5.0)  GO  TO  900 

IF  (VS.LT.5.0)  GO  TO  250 

C1=2.0*EXP(-VS) 

GO  TO  260 
250  C1=1.0/COSH(VS) 
260  CONTINUE 

IF  (VC.LT.5.0)  GO  TO  270 

C2=2.0*EXP(-VC) 

GO  TO  280 
270  C2=1.0/COSH(VC) 
280  CONTINUE 

IF  (VM.LT.5.0)  GO  TO  320 

C3=0.5*EXP(VM) 

GO  TO  330 
320  C3=C0SH(VM) 
330  CONTINUE 

FUNZ=GFUNC*TERMC*C1*C2*C3/GAMMA 

RETURN 

800    FUNZ=GFUNC*TERMC*2 . 0*EXP (GAMMA+Z) /GAMMA 
RETURN 

900    FUNZ=GFUNC*TERMC*EXP (VS+GAMMA*Z) / (COSH (VS) * GAMMA) 
RETURN 

1000  FUNZ=(  Kl/  K3)*  L3+(  Kl/  K2)*(  L1+  L2+Z) 
RETURN 
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C  BOTTOM  LAYER  CALCULATION 

C        THIS  IS  THE  BOTTOM  LAYER  CALCULATION  WHICH  IS  USED  IF  Z  FALLS 
C        INTO  THE  BOTTOM  LAYER 

C********************************************^ 
1500  IF  (GAMMA. EQ. 0.0)  GO  TO  2000 

IF  (VS.GT.5.0.AND.VC.GT.5.0.AND.VI.GT.5.0.AND. 
1  VB.GT.5.0)  GO  TO  1900 

IF(VB.GT.5.0.AND.VS.GT.5.0.AND.VC.LT.5.0.AND.VI.LT.5.0)GO  TO  2100 

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.GT.5.0.AND.VI.LT.5.0)GO  TO  2200 

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.LT.5.0.AND.VI.GT.5.0)GO  TO  2300 

IF(VB.GT.5.0.AND.VS.GT.5.0.AND.VC.GT.5.0.AND.VI.LT.5.0)GO  TO  2400 

IF(VB.GT.5.0.AND.VS.GT.5.0.AND.VC.LT.5.0.AND.VI.GT.5.0)GO  TO  2500 

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.GT.5.0.AND.VI.GT.5.0)GO  TO  2600 

IF  (VS.LT.5.0)  GO  TO  1550 

C1=2.0*EXP(-VS) 

GO  TO  1560 
1550  C1=1.0/COSH(VS) 
1560  CONTINUE 

IF  (VC.LT.5.0)  GO  TO  1570 

C2=2.0*EXP(-VC) 

GO  TO  1580 
1570  C2=1.0/COSH(VC) 
1580  CONTINUE 

IF  (VI.LT.5.0)  GO  TO  1590 

C3=2.0*EXP(-VI) 

GO  TO  1600 
1590  C3=l.O/C0SH(VI) 
1600  CONTINUE 

C4=SINH(VB) 

FUNZ=GFUNC*C1*C2*C3*C4/GAMMA 
RETURN 

1900  FUNZ=GFUNC*4 . 0*EXP (GAMMA*Z) /GAMMA 
RETURN 

2000  FUNZ=(  Kl/  K3)*(  L1+  L2+  L3+Z) 
RETURN 

2100  FUNZ=GFUNC*EXP (VB-VS) / (GAMMA*COSH (VC) *COSH (VI)  ) 
RETURN 

2200  FUNZ=GFUNC*EXP (VB-VC) / (GAMMA+COSH  (VS) *COSH (VI) ) 
RETURN 

2300  FUNZ=GFUNC*EXP (VB-VI) / (GAMMA*COSH (VS) *COSH (VC) ) 
RETURN 

2400  FUNZ=GFUNC*2 . 0*EXP (VB-VS-VC) / (GAMMA*COSH (VI) ) 
RETURN 

2500  FUNZ=GFUNC*2 . 0*EXP (VB-VS-VI) / (GAMMA*COSH (VC) ) 
RETURN 

2600  FUNZ=GFUNC*2 . 0*EXP (VB-VC-VI) / (GAMMA*COSH (VS) ) 
RETURN 
END 
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data  files  are  presented  and  discussed  to  show  both  the  increased  flexibility  of  the  input  data  and  the  actual  use  of  the 
updated  code.  Running  the  TXYZ20  code  for  one  of  the  input  files  provides  a  benchmark  for  several  machines.  The  user 
may  wish  to  run  this  example  for  the  purpose  of  comparing  the  CPU  times  involved.  The  appendix  contains  an  annotated, 
internally  documented  listing  of  the  FORTRAN  source  code  for  TXYZ20.  The  FORTRAN  source  code  (total  of  about  21 
kbytes)  and  sample  input  and  output  data  files  are  available  in  ASCII  format  using  a  number  of  transfer  vehicles.  These 
include:  standard  8  track  magnetic  tape  (ASCII,  density  =  1600,  record  =  80,  block  =  1600),  5.25-in.  (360-kbyte  and 
1.2-Mbyte)  DOS-formatted  floppy  disks,  and  electronic  mail  over  the  Internet.  The  sample  input  and  output  data  files  are 
included  so  that  the  user  can  check  the  program  for  proper  operation  as  well  as  to  become  acquainted  with  the  setup  and 
use  of  the  code.  Users  of  the  TXYZ  code  will  find  the  updated  TXYZ20  code  easy  to  use  and  should  benefit  from  the 
more  flexible  input  and  the  more  general  treatment  of  heat  sources  and  heat  sinks. 
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